

load(file="CCES_Panel_dichotomous.RData")

data2010 <- data[,grep("2010",colnames(data))]
dim(data2010)

data2012 <- data[,grep("2012",colnames(data))]
dim(data2012)

data2014 <- data[,grep("2014",colnames(data))]
dim(data2014)

###
### IRT estimator
###
source("../IRT_mixing/mixture_irt.R")

res2010 <- em_mix_irt(as.matrix(data2010), iter=30)
res2012 <- em_mix_irt(as.matrix(data2012), iter=30)
res2014 <- em_mix_irt(as.matrix(data2014), iter=30)

#save("res2010","res2012","res2014",file="cces_panel_res.RData")
#load(file="cces_panel_res.RData")

# name item parameters for items
names(res2010$irt$b) <- names(data2010)
names(res2012$irt$b) <- names(data2012)
names(res2014$irt$b) <- names(data2014)

###
### standardize the parameters, x is mean 0 SD 1
###

sdx <- sd(res2010$irt$x)
meanx <- mean(res2010$irt$x)
res2010$irt$b <- res2010$irt$b*sdx
res2010$irt$a <- res2010$irt$a +
       res2010$irt$b*meanx
res2010$irt$x <- (res2010$irt$x-meanx)/sdx

sdx <- sd(res2012$irt$x)
meanx <- mean(res2012$irt$x)
res2012$irt$b <- res2012$irt$b*sdx
res2012$irt$a <- res2012$irt$a +
       res2012$irt$b*meanx
res2012$irt$x <- (res2012$irt$x-meanx)/sdx

sdx <- sd(res2014$irt$x)
meanx <- mean(res2014$irt$x)
res2014$irt$b <- res2014$irt$b*sdx
res2014$irt$a <- res2014$irt$a +
       res2014$irt$b*meanx
res2014$irt$x <- (res2014$irt$x-meanx)/sdx

###
###
###

table(res2010$w[,1] > .5, res2012$w[,1] > .5)
table(res2012$w[,1] > .5, res2014$w[,1] > .5)

tab <- table((res2010$w[,1] > .5),(res2012$w[,1] > .5))
round(tab/sum(tab),3)

tab <- table((res2012$w[,1] > .5),(res2014$w[,1] > .5))
round(tab/sum(tab),3)


tab <- table((res2014$w[,1] > .5) == (res2012$w[,1] > .5))
tab/sum(tab)
cor((res2014$w[,1] > .5),(res2012$w[,1] > .5))
cor(res2014$w[,1],res2012$w[,1])

tab <- table((res2014$w[,1] > .5) == (res2012$w[,1] > .5) &
               (res2012$w[,1] > .5) == (res2010$w[,1] > .5) )
tab/sum(tab)

ind <- ! (res2014$w[,1] > .5 & res2012$w[,1] > .5)
cor(res2014$irt$x[ind],res2012$irt$x[ind])
ind <- res2014$w[,1] > .5 & res2012$w[,1] > .5
cor(res2014$irt$x[ind], res2012$irt$x[ind])

ind <- ! (res2010$w[,1] > .5 & res2012$w[,1] > .5)
cor(res2010$irt$x[ind],res2012$irt$x[ind])
ind <- res2010$w[,1] > .5 & res2012$w[,1] > .5
cor(res2010$irt$x[ind], res2012$irt$x[ind])


pdf("stability_of_ideal_points.pdf")

ind <- res2012$w[,1] < .5 & res2014$w[,1] < .5
plot(res2014$irt$x[ind] ~ res2012$irt$x[ind],main="",
     xlab="Ideal Point in 2012",ylab="Ideal Point in 2014",axes=F)
axis(1)
axis(2)
abline(lm(res2014$irt$x[ind] ~ res2012$irt$x[ind]))
cor(res2014$irt$x[ind],res2012$irt$x[ind])

ind <- res2012$w[,1] > .5 & res2014$w[,1] > .5
plot(res2014$irt$x[ind] ~ res2012$irt$x[ind],main="",
     xlab="Ideal Point in 2012",ylab="Ideal Point in 2014",axes=F)
axis(1)
axis(2)
abline(lm(res2014$irt$x[ind] ~ res2012$irt$x[ind]))
cor(res2014$irt$x[ind],res2012$irt$x[ind])

dev.off()

ind <- res2012$w[,1] < .5 & res2010$w[,1] < .5
cor(res2010$irt$x[ind],res2012$irt$x[ind])
cor(res2010$irt$x[! ind],res2012$irt$x[! ind])
ind <- res2012$w[,1] > .5 & res2010$w[,1] > .5
cor(res2010$irt$x[ind],res2012$irt$x[ind])


